Computing the Complete Gravitational Wavetrain from Relativistic Binary Inspiral 
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We present a new method for generating the nonhnear gravitational wavetrain horn the late in- 
spiral (pre-coalescence) phase of a binary neutron star system by means of a numerical evolution 
calculation in full general relativity. In a prototype calculation, we produce 214 wave cycles from 
corotating polytropes, representing the final part of the inspiral phase prior to reaching the ISCO. 
Our method is based on the inequality that the orbital decay timescale due to gravitational radiation 
is much longer than an orbital period and the approximation that gravitational radiation has little 
effect on the structure of the stars. We employ quasi-equilibrium sequences of binaries in circular 
orbit for the matter source in our field evolution code. We compute the gravity-wave energy flux, 
and, from this, the inspiral rate, at a discrete set of binary separations. From these data, we con- 
struct the gravitational waveform as a continuous wavetrain. Finally, we discuss the limitations of 
our current calculation, planned improvements, and potential applications of our method to other 
inspiral scenarios. 
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The inspiral and coalescence of binary neutron stars 
are among the most promising sources of gravitational 
radiation for the laser interferometer gravitational wave 
observatories LIGO, VIRGO, GEO, and TAMA, which 
should become operational in a few years. It is estimated 
§ that LIGO/VIRGO will be able to detect the last 15 
minutes of the inspiral (16,000 orbits) of a typical binary 
neutron star system. Detection and analysis of the signal 
will depend on matched filter techniques, which require 
theoretical templates of the waveforms. 

The evolution of binary neutron stars proceeds in dif- 
ferent stages. By far the longest is the initial quasi- 
equilibrium inspiral phase, during which the stars move 
in nearly circular orbits, while the separation between 
the stars slowly decreases as energy is carried away by 
gravitational radiation. The quasi-circular orbits become 
unstable at the innermost stable circular orbit (ISCO), 
where the inspiral enters a plunge and merger phase. The 
merger and coalescence of the stars happens on a dy- 
namical timescale, and produces either a black hole or a 
larger neutron star, which may collapse to a black hole 
at a later time. The final stage of the evolution is the 
ringdown phase, during which the merged object settles 
down to equilibrium. 

Two distinct approaches have commonly been em- 
ployed in general relativity to study the inspiral phase. 
Much progress has been made in post-Newtonian (PN) 
studies of compact binaries (see, e.g. g], and references 
therein). Most of these approaches, however, approxi- 
mate the stars as point sources, which neglects important 
finite-size effects for neutron stars. Also, PN expansions 
may not converge sufficiently rapidly in the strong-field 
region near the ISCO. In an independent approach, rela- 
tivistic quasi-equilibrium models of the inspiral phase of 
neutron star binaries have been constructed numerically 



(see, e.g., [|| for corotating binaries and Q for the more 
realistic case of irrotational binaries) . Binaries have also 
been evolved using fully relativistic hydrodynamics [g,|| , 
but computational constraints limit these evolutions to a 
few orbits. 

In this paper, we propose a "hydro without hydro" 
method for using the "snapshots" generated by quasi- 
equilibrium codes to create a complete wavetrain [0,p|. 
We use these snapshots to find the gravitational wave- 
form at a given separation and construct an effective ra- 
diation reaction. From this, the inspiral rate is known, 
and the entire wavetrain can be constructed. We de- 
scribe this method below. We then present the results of 
a prototype calculation, and discuss its implications. 

The fundamental assumption in the quasi-equilibrium 
approximation is that the orbital timescale is much 
shorter than the gravitational radiation reaction time- 
scale. The binary inspiral then proceeds along a sequence 
of quasi-equilibrium configurations in nearly circular or- 
bits with constant rest mass Afo- This assumption sig- 
nificantly simplifies the problem, since the hydrodynam- 
ical equations can be integrated to yield a relativistic 
Bernoulli equation |g|]. A similar approximation is rou- 
tinely adopted in stellar evolution calculations. There 
too, the evolutionary timescale is much longer than the 
hydrodynamical timescale, so that the star can safely 
be assumed to be in quasi-equilibrium on a dynamical 
timescale. 

Constructing quasi-equilibrium binary configurations 
yields their total mass-energy M{r) and their orbital an- 
gular frequency r^orb('') as a function of separation r. The 
ISCO of an evolutionary Mq = const sequence, where 
the assumption of quasi-equilibrium breaks down, can be 
identified by locating a turning point along an M versus 
r curve (see, e.g., the cross in Fig. 1). 



Given a relativistic quasi-equilibrium configuration, we 
can compute the emitted gravitational radiation by em- 
ploying the binary data for the matter source terms in 
Einstein's field equations. In particular, the gravitational 
fields can be evolved in the presence of these predeter- 
mined matter sources with a relativistic evolution code 
("hydro without hydro", M). There is no need to re- 
solve any of the hydrodynamic equations for the source, 
as we are assuming that the stars move in circular orbits 
and are little affected by the waves. We thus simply ro- 
tate the source terms on our numerical grid. Repeating 
this calculation for members of an evolutionary sequence 
at discrete separations r, we obtain the wave amplitude 
A{r) and the orbit-averaged gravitational wave luminos- 
ity < dM/dt > (r) at these separations. A fitting func- 
tion can then be used to interpolate to any intermediate 
separation (see Fig. 1). 

The inspiral is determined by computing the effective 
radiation reaction due to the radiation of mass-energy by 
gravitational waves. The energy loss < dM/dt > forces 
the binary to a closer separation, as it "slides down" the 
quasi-equilibrium M vs r curve. The inspiral rate is thus 

dr < dM/dt > 

'di ~ (dM/dr\^ ' ^ ' 

Solving this equation determines r as a function of time i, 
which then allows us to express A{t) and f7oib(i) in terms 
of t. The complete quasi-equilibrium wavetrain can then 
be constructed from 



TsKt) ^ A(t) sin( / ^G^{t')dt!) , 
Jo 



(2) 



where we expect i^GW ~ 2iloib due to the m — 2 domi- 
nance of the radiation, and where r^ is the distance be- 
tween the source and the distant observer. The validity 
of our method and its ability to produce reliable wave- 
trains is most easily demonstrated in relativistic scalar 
gravitation, where the quasi-equilibrium wavetrain can 
be compared to an exact numerical solution g] . 

In the following we present a prototype calculation 
which employs the quasi-equilibrium models of Bauni- 
garte tt.al. 0. These models assume maximal slicing 
K = A''j = 0, where Kij is the extrinsic curvature, and 
spatial conformal flatness, so that the spatial metric 7^ 
can be written as a product of the conformal factor f/; and 
a flat metric fij , 7^ = ip'^fij ■ The latter is believed to ap- 
proximately minimize the gravitational wave content of 
the spacetimes [|0| while introducing only small devia- 
tions from correct evolutionary sequences [^ . The stars 
are modeled as corotating polytropcs P = np^ '" with 
polytropic index n = \. We non-dimcnsionalizc our re- 
sults by setting k — 1. In addition to the matter profile, 
these models provide the conformal factor -0, the lapse 
function a and the shift vector p\ 

We insert these quasi-equilibrium data into our rela- 
tivistic evolution code, which is based on a conformal 
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FIG. 1. Total mass M and gravitational wave luminosity 
dM/dt versus separation for our evolutionary sequence. Both 
quantities are measured by a distant observer. The separa- 
tion is parameterized by za (see text). The points represent 
individual numerical models, and the solid lines are fits. The 
turning point dM/dzA = in the top panel locates the ISCO 
(marked by the cross). 



decomposition of the ADM equation [Q . This decompo- 
sition conveniently splits the fundamental variables into 
longitudinal quantities (1/' and K) and "wave-like" trans- 
verse quantities (the conformally related metric 7ij and 
the tracefree part of the extrinsic curvature Aij ) . We also 
introduce "conformal connection functions" F* as inde- 
pendent auxiliary quantities, which improves numerical 
stability dramatically (cf. |13[). 

We insert "by hand" the matter variables, the gauge 
variables (a and /3') and the longitudinal quantities (?/; 
and K = Q) for all times, appropriately rotated by 
dcjy/dt = fioib Pj]- The transverse variables (7^^, Aij 
and F*) are evolved dynamically, adopting the quasi- 
equilibrium data Aij , 7^^ = /^ and F* = as initial data. 
We found that we could reduce the numerical errors by 
explicitly splitting quasi-equilibrium "background" terms 
from dynamically evolved terms in the evolution equation 
for the V\ 

We evolve on a 120x120x60 zone grid, imposing equa- 
torial symmetry across the z-axis. Each star is covered by 
roughly 8 points across a diameter. We impose approx- 
imate outgoing wave boundary conditions on the trans- 
verse variables at a separation of 0.1 Agw - 0.3Agw, where 
Agw is a gravitational wavelength. We expect our outer 
boundary location to be one of the primary sources of 
error (see discussion below, and compare with |g] where 
the boundaries are similarly placed) . We monitor the vi- 
olation in the Hamiltonian and momentum constraints 
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FIG. 2. Three representative binary configurations with 
their corresponding waveforms. We show the binaries with 
greatest and smallest separation in our sequence, and one in- 
termediate member. At right, the two polarizations of gravity 
waves {h+ and hx ) along the z-axis are shown. On the z-axis, 
the two polarizations only differ by a phase. 



and find that they are aU satisfied quite well and that 
the degrees of violation do not grow significantly with 
time. We doubled the resolution for one configuration, 
and found that A changed by only « 2%. 

To extract gravitational waves, we match to even- 
parity Moncrief variables ^im (or Zerilli functions) at the 
outer part of the grid (see p5|-p^] for details). For sim- 
plicity, we focus on / = 2, m = ±2 modes, which we 
expect to be dominant. We have checked that the con- 
tribution to the energy flux of the I ^ A, ra — ±4 mode 
is very small (~ 4%), and that of the I — 2, m — Q 
mode negligible (< 10~^). Assuming I — 2, m = ±2 
dominance, we have 



^-3^«**-)^ + <**-^fl 



A 



-U 

n V 



-V'2-2 



(3) 
(4) 



where V'2±2 denotes the amplitude of V'2±2 ■ We approx- 
imate the asymptotic value of "0;^ by its value at the 
edge of our numerical grid. Due to our neglect of self- 
consistent transverse fields in the initial data, we find a 
transient noise in the waveforms during the first light- 
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FIG. 3. The amplitude A and angular frequency ficw of 
the gravitational wave signal as a function of retarded time 
t — Ts (solid line). The dotted line is the result for New- 
tonian point masses in the quadrupole approximation. The 
top label provides times in seconds for a binary of total rest 
mass Mo = 2 x I.SMq. Our sequence ends at a frequency 
Q.c,'^ = 2.36 kHz for such a binary. 



crossing time. After that time we observe the expected 
periodic waveforms originating from the periodic binary 
orbit. 

We constuct 14 different configurations along a se- 
quence of constant total rest mass Mo = 2 x 0.1 in our 
nondimensional units. The rest mass of the individual 
stars is about 55% of the maximum allowed rest mass 
of an isolated, nonrotating star. We parameterize the 
binary separation by the dimensionless quantity za, the 
ratio of the separation between the innermost points and 
the outermost points on the stars (see g). In terms 
of Za, touching binaries correspond to z^ = and in- 
finitely separated configurations to z^ = 1. In the top 
panel of Fig. 1, we plot the ADM mass M as a func- 
tion of separation, and in the bottom panel the gravita- 
tional wave luminosity dM/dt as inferred from our "hy- 
dro without hydro" integrations. Since we expect our 
quasi-equilibrium method to break down at the ISCO, 
we construct waveforms for configurations with za > 0.1, 
which is close to but safely outside the ISCO. In Fig. 2 
we illustrate three of these configurations together with 
their extracted waveforms. In Fig. 3 we show the wave 
amplitude A and the gravitational wave frequency I^gw 
as a function of retarded time. The solid line is our nu- 
merical result, which we compare with the wavetrain of 
a Newtonian point mass system as computed from the 
quadrupole approximation (dotted line). The difference 
is caused by differences in the binding energy and the 
gravitational wave luminosity due to relativity, as well as 
the numerical limitations of our calculation (see discus- 
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FIG. 4. The full -wavetrain /i+ or hx on the 2-axis as a 
function of retarded time and cycle number. The top two 
panels show enlargements of the first four and last ten cycles. 
We also indicate the gravitational strain h for a binary of total 
rest mass Mq = 2 x 1.5Mq at a separation of 100 Mpc. 



sion below). Finally, we construct the continuous gravi- 
tational wavetrain from A and Oqwi shown in Fig. 4. 

In summary, we propose a new technique to construct 
gravitational wavetrains from inspiraling relativistic bi- 
naries and demonstrate the feasability of this technique 
by a numerical example. Our scheme is extremely effi- 
cient computationally. For the computational cost of 14 
orbits, we produced the complete wavetrain covered by 
our sequence (214 cycles). There is no reason why one 
could not, with more points along the sequence, produce 
many more cycles of the inspiral. 

It is clear, however, that several improvements must 
be implemented to produce more accurate and reliable 
wavetrains. Our outer boundary conditions are crude 
and imposed at too small a separation. To gauge the 
error of our simplistic wave extraction, we measured V'22 
at 0.23Agw and at 0.35Agw and found a difference of 
roughly 10%. Gravitational waves should be extracted in 
the wave zone r^ > Agw, or more sophisticated bound- 
ary conditions and wave extraction should be applied 
(e.g. ||l^). Furthermore, the assumption of corotation is 
not realistic; using irrotational sequences is more physi- 
cal [^ . Finally, one could consider iterating this method 
by adding the radial velocity to the source and incorpo- 
rating the calculated transverse components of the field, 
as well as deviations from conformal flatness, and then 
recomputing the waves and inspiral. 

The binary black hole inspiral problem could be ap- 
proached in a similar way ||20|| . This is a prospect we are 
currently investigating. 
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